%raw_rimd_plot

stackplot = false;

% pathname = 'C:\Users\crozi\OneDrive - University Of Oregon\Documents\volcano_seismic_kilauea\Creating_txt_data\newtxtdata\';
% % tnoisestart = 0;
% % tnoiseend = 60*18;
% % teventstart = 60*22;
% % teventend = 60*22 + 60*7;
% eventpick = 'max';
% rel_tnoisestart = -60*10;
% rel_tnoiseend = -60*2;
% rel_teventstart = 60*0.25;
% rel_teventend = 60*8.25;
% eventnames = {'collapse_2018-05-17T04-00-00'};%,...
% %         'collapse_2018-05-17T13-50-00',...
% %         'collapse_2018-05-19T09-43-00',...
% %         'collapse_2018-05-20T01-43-00',...
% %         'collapse_2018-05-20T21-35-00',...
% %         'collapse_2018-05-21T10-40-00',...
% %         'collapse_2018-05-22T03-36-00',...
% %         'collapse_2018-05-23T07-43-00',...
% %         'collapse_2018-05-24T04-29-00'};

% eventnames = {'raw_explosion_2018-5-17_4-10',...
%         'raw_explosion_2018-5-17_14-0',...
%         'raw_explosion_2018-5-19_9-53',...
%         'raw_explosion_2018-5-20_1-53',...
%         'raw_explosion_2018-5-20_21-45',...
%         'raw_explosion_2018-5-21_10-50',...
%         'raw_explosion_2018-5-22_3-46',...
%         'raw_explosion_2018-5-23_7-53',...
%         'raw_explosion_2018-5-24_4-39'};

% eventnames = {'rockfall_2013-08-24T07-38-00',...
%     'rockfall_2014-07-23T20-01-00',...
%     'rockfall_2014-10-19T21-40-00',...
%     'rockfall_2015-04-25T11-50-00',...
%     'rockfall_2015-04-28T20-10-00',...
%     'rockfall_2015-05-03T23-10-00',...
%     'rockfall_2016-01-03T00-07-00',...
%     'rockfall_2016-01-08T13-40-00',...
%     'rockfall_2016-04-19T09-50-00',...
%     'rockfall_2016-08-07T07-52-00',...
%     'rockfall_2016-09-20T07-22-00',...
%     'rockfall_2016-10-10T22-15-00',...
%     'rockfall_2016-10-19T17-35-00',...
%     'rockfall_2016-11-28T21-49-00',...
%     'rockfall_2016-12-02T16-48-00'};

pathname = 'C:\Users\crozi\OneDrive - University Of Oregon\Documents\craco_waves\data\';
%tnoise = 60*14;
eventpick = 'max';
rel_tnoisestart = -60*16;
rel_tnoiseend = -60*2;
rel_teventstart = -60*2;
rel_teventend = 60*12;
eventnames = {'rockfall_2013-08-24T07-33-00'};%,...
%     'rockfall_2014-07-23T19-56-00',...
%     'rockfall_2014-10-19T21-35-00',...
%     'rockfall_2015-04-25T11-45-00',...
%     'rockfall_2015-04-28T20-05-00',...
%     'rockfall_2015-05-03T23-05-00',...
%     'rockfall_2016-01-03T00-02-00',...
%     'rockfall_2016-01-08T13-35-00',...
%     'rockfall_2016-04-19T09-45-00',...
%     'rockfall_2016-08-07T07-47-00',...
%     'rockfall_2016-09-20T07-17-00',...
%     'rockfall_2016-10-10T22-10-00',...
%     'rockfall_2016-10-19T17-30-00',...
%     'rockfall_2016-11-28T21-44-00',...
%     'rockfall_2016-12-02T16-43-00'};

okadacell = struct('eventname',eventnames);

datestrings  = cell(size(eventnames));
for i = 1:length(datestrings)
    datestrings{i} = eventnames{i}(10:end);
end
datetimes = datetime(datestrings,'InputFormat','yyyy-MM-dd''T''HH-mm-ss');
    

for k = 1:length(okadacell)
    %%%%import seismic data
    %import HVO station locations
    stations = {'KKO'; 'BYL'; 'HAT'; 'SBL'; 'UWB'; 'NPT'; 'OBL'; 'UWE'; 'SDH'; 'WRM'; 'RIMD'; 'NPB'; 'SRM'; 'PUHI'};
    %stations = {'KKO'; 'BYL'; 'HAT'; 'SBL'; 'UWB'; 'NPT'; 'OBL'; 'UWE'; 'SDH'; 'WRM';'RIMD'};
    stationindex = true(size(stations));
    for j = 1:length(stations)
        if ~exist(strcat(pathname,'halemaumau_',eventnames{k},'_',stations{j},'.txt'),'file')
            stationindex(j) = false;
        end
    end
    stations = stations(stationindex);
    locations = cell(size(stationindex));
    locations{1} = struct('station','KKO', 'latlon', [19.39781 -155.26600 1146.0], 'misfit_weight',1);
    locations{2} = struct('station','BYL', 'latlon', [19.41200 -155.26000 1059.0],'misfit_weight',0);
    locations{3} = struct('station','HAT', 'latlon', [19.42300 -155.26100 1082.0],'misfit_weight',0);
    locations{4} = struct('station','SBL', 'latlon', [19.42700 -155.26800 1084.0],'misfit_weight',1);
    locations{5} = struct('station','UWB', 'latlon', [19.42469 -155.27800 1091.0],'misfit_weight',1);
    locations{6} = struct('station','NPT', 'latlon', [19.41200 -155.28100 1115.0],'misfit_weight',1);
    locations{7} = struct('station','OBL', 'latlon', [19.41750 -155.28400 1107.0],'misfit_weight',1);
    locations{8} = struct('station','UWE', 'latlon', [19.42111 -155.29300 1240.0],'misfit_weight',1);
    locations{9} = struct('station','SDH', 'latlon', [19.38950 -155.29400 1133.0],'misfit_weight',1);
    locations{10} = struct('station','WRM', 'latlon', [19.40650 -155.30000 1163.0],'misfit_weight',1);
    locations{11} = struct('station','RIMD', 'latlon', [19.39531 -155.27400 1128.0],'misfit_weight',1);
    locations{12} = struct('station','NPB', 'latlon', [0 0 0],'misfit_weight',1);
    locations{13} = struct('station','SRM', 'latlon', [0 0 0],'misfit_weight',1);
    locations{14} = struct('station','PUHI', 'latlon', [19.385510 -155.251310 1079],'misfit_weight',1);
    %convert lat-lon locations to UTM E-N coordinates
    zone = utmzone(locations{1}.latlon(1),locations{1}.latlon(2));
    utmstruct = defaultm('utm');
    utmstruct.zone = zone;
    utmstruct.geoid = wgs84Ellipsoid;
    utmstruct = defaultm(utmstruct);
    for i = 1:length(locations)
        locations{i}.loc = [0 0 0];
       [locations{i}.loc(1), locations{i}.loc(2), locations{i}.loc(3)] = ...
           mfwdtran(utmstruct,locations{i}.latlon(1),locations{i}.latlon(2),locations{i}.latlon(3)); 
    end

    %import waveforms and instrument responses
    okadacell(k).s = cell(size(stations));
    count = 0;
    for i = 1:length(locations)
        if ismember(locations{i}.station,stations)
            count = count + 1;
            okadacell(k).s{count} = struct('station',locations{i}.station);
            okadacell(k).s{count}.misfit_weight = locations{i}.misfit_weight;
            %waveforms
            data = csvread(strcat(pathname,'halemaumau_',eventnames{k},'_',okadacell(k).s{count}.station,'.txt'),1);
            okadacell(k).s{count}.time_raw = data(:,1);
            okadacell(k).s{count}.east_raw = data(:,2);
            okadacell(k).s{count}.north_raw = data(:,3);
            okadacell(k).s{count}.vert_raw = data(:,4);
            if mod(length(okadacell(k).s{count}.time_raw),2) == 1
                okadacell(k).s{count}.time_raw = okadacell(k).s{count}.time_raw(1:end-1);
                okadacell(k).s{count}.east_raw = okadacell(k).s{count}.east_raw(1:end-1);
                okadacell(k).s{count}.north_raw = okadacell(k).s{count}.north_raw(1:end-1);
                okadacell(k).s{count}.vert_raw = okadacell(k).s{count}.vert_raw(1:end-1);
            end
            okadacell(k).s{count}.loc_east = locations{i}.loc(1);
            okadacell(k).s{count}.loc_north = locations{i}.loc(2);
            okadacell(k).s{count}.loc_vert = locations{i}.loc(3);
%             %instrument responses
%             data = csvread(strcat(s{count}.station,'_vel_response.txt'),1);
%             s{count}.raw_response = data(:,2)+1i*data(:,3);
%             s{count}.raw_f_response = data(:,1);
        end
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%plotting
window = 120;%seconds
toverlap = 5; %seconds
maxfreqspec = 10;
maxfreqpeaks = 0.15;
for i = 1:length(datetimes)
    xticklist = 0:60:6000;
    xticklabellist = datestr(datetimes(i) + seconds(xticklist),'HH:MM');
    xticklabellist(1:2:end,:) = ' ';
    
  
    figure('name','separate')
    
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %waveforms
    
    subplot(3,1,1)
    for k = 1:length(okadacell(i).s)
        if strcmp(okadacell(i).s{k}.station, 'RIMD')
            index = k;
        end
    end
    plot(okadacell(i).s{index}.time_raw,...
        okadacell(i).s{index}.vert_raw - mean(okadacell(i).s{index}.vert_raw), 'k')
    hold on
    grid on
    title(strcat('Rockfall',{' '}, datestrings{i}(1:10)))
    xticks(xticklist)
    xtickangle(0)%, grid minor
    xticklabels(xticklabellist)
    ylabel('Counts')
    xlabel('Time (UTC)')
%     xlim([tevent,max(okadacell(i).s{index}.time_raw)-floor(window/2)])
    ylim([min(okadacell(i).s{index}.vert_raw - mean(okadacell(i).s{index}.vert_raw)),...
        max(okadacell(i).s{index}.vert_raw - mean(okadacell(i).s{index}.vert_raw))])
    
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %bandpass filter
    Fstop1 = 1/120;    % First Stopband Frequency
    Fpass1 = 1/90;    % First Passband Frequency
    Fpass2 = 1/8;     % Second Passband Frequency
    Fstop2 = 1/4;     % Second Stopband Frequency
    Dstop1 = 10;     % First Stopband Attenuation
    Dpass  = 2;  % Passband Ripple
    Dstop2 = 10;     % Second Stopband Attenuation
    flag   = 'scale';           % Sampling Flag

    fsraw = 1/(okadacell(i).s{index}.time_raw(2)...
        - okadacell(i).s{index}.time_raw(1));
    
    filt = designfilt('bandpassiir','Samplerate',fsraw,...
        'StopbandFrequency1',Fstop1,...
        'PassbandFrequency1',Fpass1,...
        'PassbandFrequency2',Fpass2,...
        'StopbandFrequency2',Fstop2,...
        'StopbandAttenuation1',Dstop1,...
        'PassbandRipple',Dpass,'StopbandAttenuation2',Dstop2,...
        'MatchExactly', 'passband','DesignMethod','butter');
    delay = 0;%floor(mean(grpdelay(filt))); %use if using filter instead of filtfilt
    
    edgetaperwavelength = 400/Fpass2;
    edgetaper = (tukeywin(length(okadacell(i).s{index}.time_raw),...
        edgetaperwavelength...
        /(okadacell(i).s{index}.time_raw(end)-okadacell(i).s{index}.time_raw(1))));

    okadacell(i).s{index}.vert_band = filtfilt(filt,...
        [(okadacell(i).s{index}.vert_raw - mean(okadacell(i).s{index}.vert_raw))...
        .*edgetaper; zeros(delay,1)]);
    okadacell(i).s{index}.vert_band = (okadacell(i).s{index}.vert_band(delay+1:end));
    plot(okadacell(i).s{index}.time_raw, okadacell(i).s{index}.vert_band, 'r')
    
    legend({'RIMD Vertical','8-90 s Band RIMD Vertical'})
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %plot spectrograms
    subplot(3,1,2)
    colormap('jet')
    Ts  = (okadacell(i).s{1}.time_raw(2)-okadacell(i).s{1}.time_raw(1));
    noverlap = floor(window/Ts)-floor(toverlap/Ts);
    for j = 1:length(okadacell(i).s)
        if j == 1
            [sspec,fspec,tspec] = spectrogram(okadacell(i).s{j}.vert_raw, floor(window/Ts),noverlap,floor(window/Ts),1/Ts);
        else
            sspec = sspec + spectrogram(okadacell(i).s{j}.vert_raw, floor(window/Ts),noverlap,floor(window/Ts),1/Ts);
        end
        sspec = sspec + spectrogram(okadacell(i).s{j}.east_raw, floor(window/Ts),noverlap,floor(window/Ts),1/Ts);
        sspec = sspec + spectrogram(okadacell(i).s{j}.north_raw, floor(window/Ts),noverlap,floor(window/Ts),1/Ts);
    end
    sspec_plot = log10(abs((sspec(fspec<maxfreqspec & fspec>fspec(2),:))));
    h = pcolor(tspec,1./fspec(fspec<maxfreqspec & fspec>fspec(2)),sspec_plot);
    set(h, 'EdgeColor', 'none');
    set(gca,'ydir','normal')
    set(gca,'YScale', 'log');
%     imagesc(t,f(f<10),log10(abs((s(f<10,:)))))
%     set(gca,'ydir','normal')
%     ylim(1./fliplr([1/window,2]))
    %ylim([6 12])
    caxis([min(sspec_plot(:)) + 1/3*(max(sspec_plot(:)) - min(sspec_plot(:))),...
        max(sspec_plot(:)) - 1/20*(max(sspec_plot(:)) - min(sspec_plot(:)))])
    title('Stacked Seismic Counts from all Stations')
    xlabel('Time (UTC)')
    ylabel('Period (s)')
%     ylabel('Frequency (Hz)')
    cbar = colorbar;
    cbar.Position = cbar.Position.*[0,1,0,1] + [0.91,0,0.01,0];
    ylabel(cbar,'Log_{10} Counts')
    xticks(0:60:6000)
    %yticks(.2:.2:2)
    %yticks([[1]*10^-2,[1]*10^-1,[1]*10^0,[1]*10^1,[1]*10^2])
    grid on
    set(gca, 'YMinorTick','on', 'YMinorGrid','on')
    set(gca,'layer','top')
    xticks(xticklist)
    xtickangle(0)%, grid minor
    xticklabels(xticklabellist)
%     xlim([tevent,max(okadacell(i).s{index}.time_raw)-floor(window/2)])
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %PSD
    subplot(3,1,3)
    %}

    Fs = 1/(okadacell(i).s{1}.time_raw(2) - okadacell(i).s{1}.time_raw(1));

    if strcmp(eventpick,'max')
        %pick based on spectrogram max
        fmax_ind = find(fspec>maxfreqspec,1);
        sspec_sum = sum(abs(sspec(3:fmax_ind,:)),1);
        [~,ievent] = max(sspec_sum);
        tevent = tspec(ievent);
        tnoisestart = tevent + rel_tnoisestart;
        tnoiseend = tevent + rel_tnoiseend;
        teventstart = tevent + rel_teventstart;
        teventend = tevent + rel_teventend;
    end
    
    inoisestart = find(okadacell(i).s{j}.time_raw > tnoisestart, 1);
    inoiseend = find(okadacell(i).s{j}.time_raw > tnoiseend, 1);
    if mod(inoiseend - inoisestart,2) == 0
        inoisestart = inoisestart + 1;
    end
    Nnoise = length(okadacell(i).s{1}.vert_raw(inoisestart:inoiseend));
    freqnoise = (0:Fs/Nnoise:Fs/2).';

    ieventstart = find(okadacell(i).s{j}.time_raw > teventstart,1);
    ieventend = find(okadacell(i).s{j}.time_raw > teventend,1);
    if mod(ieventstart - ieventend,2) == 0
        ieventstart = ieventstart + 1;
    end
    Nevent = length(okadacell(i).s{j}.vert_raw(ieventstart:ieventend));
    freqevent = (0:Fs/Nevent:Fs/2).';

    taperwavelength = 20;
    tapernoise = tukeywin(Nnoise, taperwavelength*Fs);
    taperevent = tukeywin(Nevent, taperwavelength*Fs);

    for j = 1:length(okadacell(i).s)        
        %noise
        xdftnoise = fft(detrend(okadacell(i).s{j}.vert_raw(inoisestart:inoiseend)).*tapernoise);
        xdftnoise = xdftnoise(1:Nnoise/2+1);
        psdnoise = (1/(Fs*Nnoise)) * abs(xdftnoise).^2;
        %psdnoise = abs(xdftnoise).^2;
        psdnoise(2:end-1) = 2*psdnoise(2:end-1);
        if j == 1       
            snoise = psdnoise;
            [swelchnoise, freqwelchnoise] = pwelch(detrend(okadacell(i).s{j}.vert_raw(inoisestart:inoiseend)).*tapernoise,[],[],[],Fs);
            [smtmnoise, freqmtmnoise] = pmtm(detrend(okadacell(i).s{j}.vert_raw(inoisestart:inoiseend)).*tapernoise,[],[],Fs);
        else
            swelchnoise = swelchnoise + pwelch(detrend(okadacell(i).s{j}.vert_raw(inoisestart:inoiseend)).*tapernoise,[],[],[],Fs);
            smtmnoise = smtmnoise + pmtm(detrend(okadacell(i).s{j}.vert_raw(inoisestart:inoiseend)).*tapernoise,[],[],Fs);
            snoise = snoise + psdnoise;
        end
        xdftnoise = fft(detrend(okadacell(i).s{j}.east_raw(inoisestart:inoiseend)).*tapernoise);
        xdftnoise = xdftnoise(1:Nnoise/2+1);
        psdnoise = (1/(Fs*Nnoise)) * abs(xdftnoise).^2;
        %psdnoise = abs(xdftnoise).^2;
        psdnoise(2:end-1) = 2*psdnoise(2:end-1);
        snoise = snoise + psdnoise;
        swelchnoise = swelchnoise + pwelch(detrend(okadacell(i).s{j}.east_raw(inoisestart:inoiseend)).*tapernoise,[],[],[],Fs);
        smtmnoise = smtmnoise + pmtm(detrend(okadacell(i).s{j}.east_raw(inoisestart:inoiseend)).*tapernoise,[],[],Fs);
        
        xdftnoise = fft(detrend(okadacell(i).s{j}.north_raw(inoisestart:inoiseend)).*tapernoise);
        xdftnoise = xdftnoise(1:Nnoise/2+1);
        psdnoise = (1/(Fs*Nnoise)) * abs(xdftnoise).^2;
        %psdnoise = abs(xdftnoise).^2;
        psdnoise(2:end-1) = 2*psdnoise(2:end-1);
        snoise = snoise + psdnoise;
        swelchnoise = swelchnoise + pwelch(detrend(okadacell(i).s{j}.north_raw(inoisestart:inoiseend)).*tapernoise,[],[],[],Fs);
        smtmnoise = smtmnoise + pmtm(detrend(okadacell(i).s{j}.north_raw(inoisestart:inoiseend)).*tapernoise,[],[],Fs);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %signal    
        xdftevent = fft(detrend(okadacell(i).s{j}.vert_raw(ieventstart:ieventend)).*taperevent);
        xdftevent = xdftevent(1:Nevent/2+1);
        psdevent = (1/(Fs*Nevent)) * abs(xdftevent).^2;
        %psd = abs(xdft).^2;
        psdevent(2:end-1) = 2*psdevent(2:end-1);
        if j == 1
            sevent = psdevent;
            [swelchevent, freqwelchevent] = pwelch(detrend(okadacell(i).s{j}.vert_raw(ieventstart:ieventend)).*taperevent,[],[],[],Fs);
            [smtmevent, freqmtmevent] = pmtm(detrend(okadacell(i).s{j}.vert_raw(ieventstart:ieventend)).*taperevent,[],[],Fs);
        else
            sevent = sevent + psdevent;
            swelchevent = swelchevent + pwelch(detrend(okadacell(i).s{j}.vert_raw(ieventstart:ieventend)).*taperevent,[],[],[],Fs);
            smtmevent = smtmevent + pmtm(detrend(okadacell(i).s{j}.vert_raw(ieventstart:ieventend)).*taperevent,[],[],Fs);
        end
        xdftevent = fft(detrend(okadacell(i).s{j}.east_raw(ieventstart:ieventend)).*taperevent);
        xdftevent = xdftevent(1:Nevent/2+1);
        psdevent = (1/(Fs*Nevent)) * abs(xdftevent).^2;
        %psd = abs(xdft).^2;
        psdevent(2:end-1) = 2*psdevent(2:end-1);
        sevent = sevent + psdevent;
        swelchevent = swelchevent + pwelch(detrend(okadacell(i).s{j}.east_raw(ieventstart:ieventend)).*taperevent,[],[],[],Fs);
        smtmevent = smtmevent + pmtm(detrend(okadacell(i).s{j}.east_raw(ieventstart:ieventend)).*taperevent,[],[],Fs);
        
        xdftevent = fft(detrend(okadacell(i).s{j}.north_raw(ieventstart:ieventend)).*taperevent);
        xdftevent = xdftevent(1:Nevent/2+1);
        psdevent = (1/(Fs*Nevent)) * abs(xdftevent).^2;
        %psd = abs(xdft).^2;
        psdevent(2:end-1) = 2*psdevent(2:end-1);
        sevent = sevent + psdevent;
        swelchevent = swelchevent + pwelch(detrend(okadacell(i).s{j}.north_raw(ieventstart:ieventend)).*taperevent,[],[],[],Fs);
        smtmevent = smtmevent + pmtm(detrend(okadacell(i).s{j}.north_raw(ieventstart:ieventend)).*taperevent,[],[],Fs);
    end
    
    loglog(1./freqnoise, snoise, 'c')
    hold on
    loglog(1./freqevent, sevent, 'b')
    loglog(1./freqwelchnoise, swelchnoise, 'm')
    loglog(1./freqwelchevent, swelchevent, 'r')
    loglog(1./freqmtmnoise, smtmnoise, 'y')
    loglog(1./freqmtmevent, smtmevent, 'g')
%     xlim(1./fliplr([min(freq),max(freq)]))
    xlim(1./fliplr([freqevent(4),10]))
    endin = find(freqevent>10);
    ylim([min(snoise(4:endin)),max(sevent(4:endin))*7])
    grid on
    title('Periodogram')
    ylabel('Power (Counts)')
    %xlabel('Frequency (Hz)')
    xlabel('Period (s)')
    legend({'Preceding Noise','Signal'},...
        'Location','southeast','Autoupdate','off')
    
    %find peaks
%     figure('name','peaks')
    %plot(freqmtm(2:end), log10(smtm(2:end)), 'k')
    ssmooth = smoothdata(log10(sevent(2:end)),1,'gaussian',3);
    maxind = find(freqevent >= maxfreqpeaks);
    snoisesmooth = smoothdata(log10(snoise(2:end)),1,'gaussian',9);
    snoisesmooth = interp1(freqnoise(2:end),snoisesmooth,freqevent(2:end),...
        'moving','extrap');
%     [pks,pkf,width,prom] = ...
%         findpeaks(ssmooth(1:maxind-1)...
%         ./snoisesmooth(1:maxind-1),...
%             freq(2:maxind),'SortStr','none',...
%         'MinPeakHeight', 1.04,...
%         'MinPeakProminence', 0.05,...
%         'MinPeakDistance', 0.005,...
%         'Annotate','extents');
    [pks,pkf,width,prom] = ...
        findpeaks(ssmooth(1:maxind-1)...
        ./snoisesmooth(1:maxind-1),...
            freqevent(2:maxind),'SortStr','none',...
        'MinPeakHeight', 1.025,...
        'MinPeakProminence', 0.02,...
        'MinPeakDistance', 0.003,...
        'Annotate','extents');
    %set(gca,'xscale','log')
    if ~isempty(pks)
        edges = [freqevent(4);pkf;maxfreqpeaks];
        midedges = (edges(2:end) + edges(1:end-1))/2;
        flist = NaN(size(pks));
        %Q1list = NaN(size(pks));
        %Q2list = NaN(size(pks));
        Qlist = NaN(size(pks));
        Wlist = NaN(size(pks));
        slist = NaN(size(pks));
        txtcell = cell(size(pks));
        for n = 1:length(pks)
            stemp = sevent;
            stemp([1,2,3]) = 0;
%             stemp(freq<=edges(n) | freq>=edges(n+2)) =...
%                 stemp(freq<=edges(n) | freq>=edges(n+2))/1E5;
            stemp(freqevent<=midedges(n) | freqevent>=midedges(n+1)) =...
                stemp(freqevent<=midedges(n) | freqevent>=midedges(n+1))/1E5;
            tempwindow = gmdistribution(pkf(n),0.005/200);
            tempwindow = pdf(tempwindow,freqevent);
            %[f1, Q1] = full_width_half_max(freq, stemp);
            %[~, Q2] = Robinson_and_Clegg(freq, stemp);
            %Q1list(n) = Q1;
            %Q2list(n) = Q2;
            [~,calcQi] = max(stemp.*tempwindow);
            calcQGmax = sevent(calcQi);
            calcQf0 = freqevent(calcQi); % resonance frequency
            calcQscale = sqrt(2)/2;
            calcQimin = find(sevent(1:calcQi  )<=calcQscale*calcQGmax,1,'last');
            calcQimax = find(sevent(calcQi:end)<=calcQscale*calcQGmax,1,'first')+calcQi-1;
            calcQfmin = freqevent(calcQimin);
            calcQfmax = freqevent(calcQimax);
            calcQ = calcQf0/(calcQfmax-calcQfmin); % quality factor

%             %debug plots
%             figure()
%             semilogy(freq,tempwindow)
%             xlim([0,1]),ylim(10.^[-10,20]),hold on
%             semilogy(freq,s)
%             semilogy(freq,stemp.*tempwindow)
%             semilogy(freq,stemp)

            [~,f3in] = max(stemp.*tempwindow);
            slist(n) = sevent(f3in);
            f3 = freqevent(f3in);
            flist(n) = f3;
            txtcell{n} = strcat('T=',num2str(1/f3,'%.2f'),'s');
            valid = false;
            if ~isempty(calcQimin) && ~isempty(calcQimax)
                Qlist(n) = calcQ;
                Wlist(n) = (calcQfmax-calcQfmin);
                txtcell{n} = strcat(txtcell{n},'  Q=',...
                    num2str(Qlist(n),'%.2f'));
%             else
%                 txtcell{n} = strcat(txtcell{n},'  Q=',...
%                     num2str(Qlist(n),'%.1f'));
            end
        end

        %plot peaks
        scatter(1./flist,slist,36,'k','*')
        ytxt = max(6*sevent(4:endin));
        h = text(1./(flist/1.1), ytxt+0*flist, txtcell, 'Color', 'k');
        set(h,'Rotation',-90);
        for n = 1:length(slist)
            plot([1/flist(n),1/flist(n)],...
                [ytxt,slist(n)],'k')
        end
    end
    
    
    %
    
    
    %update x axis
    subplot(3,1,1)
    xlim([teventstart,teventend])
    subplot(3,1,2)
    xlim([teventstart,teventend])
end


%{
if stackplot

    %cross correlate RIMD
    rshifts = zeros(size(datetimes));
    summed_RIMDnorth = full_RIMDnorths{1};
    for i = 2:length(datetimes)
        [acor,lag] = xcorr(full_RIMDnorths{1},full_RIMDnorths{i});
        [~,in] = max(abs(acor));
        rshifts(i) = lag(in);

        if rshifts(i) > 0
            full_RIMDnorths{i} = [zeros(rshifts(i),1);full_RIMDnorths{i}(1:end-rshifts(i))];
            full_RIMDverts{i} = [zeros(rshifts(i),1);full_RIMDverts{i}(1:end-rshifts(i))];
                okadacell(i).vert_raw = [zeros(rshifts(i),1);okadacell(i).vert_raw(1:end-rshifts(i))];
                okadacell(i).east_raw = [zeros(rshifts(i),1);okadacell(i).east_raw(1:end-rshifts(i))];
                okadacell(i).north_raw = [zeros(rshifts(i),1);okadacell(i).north_raw(1:end-rshifts(i))];
        elseif rshifts(i) < 0
            full_RIMDnorths{i} = [full_RIMDnorths{i}(1-rshifts(i):end);zeros(-rshifts(i),1)];
            full_RIMDverts{i} = [full_RIMDverts{i}(1-rshifts(i):end);zeros(-rshifts(i),1)];
                okadacell(i).vert_raw = [okadacell(i).vert_raw(1-rshifts(i):end);zeros(-rshifts(i),1)];
                okadacell(i).east_raw = [okadacell(i).east_raw(1-rshifts(i):end);zeros(-rshifts(i),1)];
                okadacell(i).north_raw = [okadacell(i).north_raw(1-rshifts(i):end);zeros(-rshifts(i),1)];
        end
        summed_RIMDnorth = summed_RIMDnorth + full_RIMDnorths{i};
    end

    figure('name','RIMDnorth separate')
    hold on
    %legendcell = {};
    yspace = max(max(abs([full_RIMDnorths{:}])))*1.5;
    for i = 1:length(full_RIMDnorths)
        %subplot(1,2*length(full_iso),[1+2*(i-1),2+2*(i-1)])
        if i > 9
            eventscale = 1;
        else
            eventscale = 1;
        end
        plot(okadacell(i).time_raw, full_RIMDnorths{i}*eventscale + yspace*((length(full_RIMDnorths)+1-i)-1),'g')
        %legendcell{i} = (okadacell(i).eventname);
        eventname = eventnames{i};
        eventname(eventname == '_') = ' ';
        %text(210,yspace/5+yspace*((length(full_RIMDnorths)+1-i)-1), strcat(eventname(14:end),'  scale x',num2str(eventscale)));
        text(210,yspace/3+yspace*((length(full_RIMDnorths)+1-i)-1), eventname(14:end))
        grid on
        if i == 1 %length(full_iso)
        xlabel('cross-correlated time (s)')
        ylabel('counts')
        end
    end
    title('RIMD North Counts')
    ylim([-yspace/2, yspace*(length(full_RIMDnorths)-0.5)])
    yticks(0:yspace:20*yspace)
    xticks(0:60:10000)
    %xlim([min(okadacell(i).s{1}.time_raw),max(okadacell(i).s{1}.time_raw)])
    xlim([200,500])

    figure('name','RIMDvert separate')
    hold on
    %legendcell = {};
    yspace = max(max(abs([full_RIMDverts{:}])))*1.5;
    for i = 1:length(full_RIMDverts)
        %subplot(1,2*length(full_iso),[1+2*(i-1),2+2*(i-1)])
        if i > 9
            eventscale = 1;
        else
            eventscale = 1;
        end
        plot(okadacell(i).time_raw, full_RIMDverts{i}*eventscale + yspace*((length(full_RIMDverts)+1-i)-1),'g')
        %legendcell{i} = (okadacell(i).eventname);
        eventname = eventnames{i};
        eventname(eventname == '_') = ' ';
        %text(210,yspace/5+yspace*((length(full_RIMDverts)+1-i)-1), strcat(eventname(14:end),'  scale x',num2str(eventscale)));
        text(210,yspace/3+yspace*((length(full_RIMDverts)+1-i)-1), eventname(14:end))
        grid on
        if i == 1 %length(full_iso)
        xlabel('cross-correlated time (s)')
        ylabel('counts')
        end
    end
    title('RIMD Vert Counts')
    ylim([-yspace/2, yspace*(length(full_RIMDverts)-0.5)])
    yticks(0:yspace:20*yspace)
    xticks(0:60:10000)
    %xlim([min(okadacell(i).s{1}.time_raw),max(okadacell(i).s{1}.time_raw)])
    xlim([200,500])

end
%}

